#include "model.h"

using namespace std;
using namespace mu;

//------------------------------------------------------------------------------
Model::Model(Mesh* mymesh,FsContainer* myfs,SaveOutputGmsh* myoutput,ClassParameter& param){

  order     = param.getInteger("Order");
  mesh      = mymesh;
  fs        = myfs;
  output    = myoutput;
  graphLap  = false;
  lumpedM   = false;
  flux      = FluxFactory::CreateNewFlux(mymesh,myfs,param);
  string gL = param.getToken("GraphLaplacian");
  if(gL=="Yes") graphLap=true;
  string lM = param.getToken("LumpedMass");
  if(lM=="Yes") lumpedM=true;
  
  // --- create containers -----------------------------------------------------
    
  int nbFQP    = (*fs)[0]->getNbFQP();
  int nbFNP    = (*fs)[0]->getNbFNP();
  int nbFaces  = mesh->getNbIntFaces();
  int nbBFaces = mesh->getNbBFaces();
  
  int eT=0;
  for(FsContainer::const_iterator fsIt = fs->begin();fsIt!=fs->end();++fsIt,eT++){
    int nbENP  = (*fsIt)->getNbENP();
    int nbEQP  = (*fsIt)->getNbEQP();
    int nbElmt = mesh->getNbElmt(eT);
    UEQP.push_back(DGMatrix(nbEQP,nbElmt));
    SolContainer  tmp; SolContainer dtmp;
    for(int dim=0;dim<2;dim++){
       tmp.push_back(DGMatrix(nbEQP,nbElmt)); dtmp.push_back(DGMatrix(nbEQP,nbElmt));
    }
    FUEQP.push_back(tmp); dFUEQP.push_back(dtmp);
    bE.push_back(DGMatrix(nbENP,nbENP));
    nuE.push_back(DGVector(nbElmt));
  }
  
  for(int dim=0;dim<2;dim++){
    UFNP.push_back(DGMatrix(nbFNP,nbFaces));
    UFQP.push_back(DGMatrix(nbFQP,nbFaces));
    UBQP.push_back(DGMatrix(nbFQP,nbBFaces));
    SolContainer tmpFN;
    SolContainer tmpF; SolContainer dtmpF;
    SolContainer tmpB; SolContainer dtmpB;
    for(int LR=0;LR<2;LR++){
      tmpFN.push_back(DGMatrix(nbFNP,nbFaces));
      tmpF.push_back(DGMatrix(nbFQP,nbFaces));  dtmpF.push_back(DGMatrix(nbFQP,nbFaces));
      tmpB.push_back(DGMatrix(nbFQP,nbBFaces)); dtmpB.push_back(DGMatrix(nbFQP,nbBFaces));
    }
    FUFNP.push_back(tmpFN);
    FUFQP.push_back(tmpF); dFUFQP.push_back(dtmpF);
    FUBQP.push_back(tmpB); dFUBQP.push_back(dtmpB);
  }
  
  bF  = DGMatrix(2*nbFNP,2*nbFNP);
  nuF = DGVector(nbFaces);
  
  setBilinearForm();
  
}


//------------------------------------------------------------------------------
double Model::getCFL(double dt){                                                // CFL = dt/h * Beta * (2p+1) 
  
  return dt*flux->getMaxBetaOverHmin()*(2*order+1);                     
  
}
//------------------------------------------------------------------------------
void Model::computeMaxPrinciple(const SolContainer* U,DGMatrix* MPminmax){
  
  int eT=0;
  for(SolContainer::const_iterator eTIt=U->begin();eTIt!=U->end();++eTIt,eT++){
    
    int nbENP  = (*fs)[eT]->getNbENP();
    int nbElmt = mesh->getNbElmt(eT);
            
    // --- Loop over the elements ---
    
    #pragma omp parallel for
    for(int eID=0;eID<nbElmt;eID++){
      
      double minENP = (*U)[eT](0,eID);
      double maxENP = (*U)[eT](0,eID);
      
      for(int i=0;i<nbENP;i++){
        if(minENP>(*U)[eT](i,eID)) minENP=(*U)[eT](i,eID);
        if(maxENP<(*U)[eT](i,eID)) maxENP=(*U)[eT](i,eID);        
      }
      
      for(int f=0;f<mesh->getNbElmtFaces(eT);f++){
        
        int fID  = mesh->getElmtFaceID(eT,eID,f);
        int LorR = mesh->getElmtFaceID(eT,eID,f,1);
        if(LorR==3) continue;                                                   //ATTENTION: BC is not yet included
        
        int R        ;LorR==0? R=1:R=0;
        int eTR      = mesh->getIntFacesTag(fID,3-3*LorR);                               
        int eIDR     = mesh->getIntFacesTag(fID,4-3*LorR); 
        int eFNumR   = mesh->getIntFacesTag(fID,5-3*LorR);
        const DGIntMatrix* indexR = (*fs)[eTR]->getFRefNodesIndex();

        for(int i=0;i<(*fs)[eTR]->getNbFNP();i++){
          if(minENP>(*U)[eTR]((*indexR)(eFNumR,i),eIDR)) minENP=(*U)[eTR]((*indexR)(eFNumR,i),eIDR);
          if(maxENP<(*U)[eTR]((*indexR)(eFNumR,i),eIDR)) maxENP=(*U)[eTR]((*indexR)(eFNumR,i),eIDR);
        }

      }

      // ---- check MP ---------
      
      if(minENP<(*MPminmax)(eID,0) || maxENP>(*MPminmax)(eID,1)){
        (*MPminmax)(eID,2)=1;
      }else{
        (*MPminmax)(eID,2)=0;
      } 
       
      (*MPminmax)(eID,0) = minENP;
      (*MPminmax)(eID,1) = maxENP;
         
    }
    
  }
  
  
}
//------------------------------------------------------------------------------
void Model::computeResidual(const SolContainer* U, SolContainer* Res){
  
  computeUQP(U);
  computeFUQP();
  
  computeVolumeTerm(Res);
  computeInterfaceTerm(Res);
  computeBoundaryTerm(Res);
  
  if(graphLap){
    computeViscosity(U);
    computeArtViscousVolumeTerm(U,Res);
    computeArtViscousInterfaceTerm(U,Res);
  }
  
  multiplyMass(Res);
  
}
//------------------------------------------------------------------------------
void Model::multiplyMass(SolContainer* Res){
  
  int eT = 0;
  for(SolContainer::const_iterator eTIt=Res->begin();eTIt!=Res->end();++eTIt,eT++){
    DGMatrix tmpres((*Res)[eT].size1(),(*Res)[eT].size2());
    if(lumpedM) {tmpres.blas_dgemm_omp_col(*(*fs)[eT]->getInvLumpedMass(),(*Res)[eT],1,0);}
    else        {tmpres.blas_dgemm_omp_col(*(*fs)[eT]->getInvMass(),(*Res)[eT],1,0);}
    (*Res)[eT].memcpy(tmpres);
  }
  
}
//------------------------------------------------------------------------------
void Model::computeUQP(const SolContainer* U){
  
  // --- get solution at elements quadrature points ----------------------------
  
  int eT = 0;
  for(SolContainer::const_iterator eTIt=U->begin();eTIt!=U->end();++eTIt,eT++){
    UEQP[eT].blas_dgemm_omp_col(*(*fs)[eT]->getPhiEQP(),*eTIt,1.,0.);                   //TO BE OPENMP
  }
  
  // --- get solution at interfaces quadrature points --------------------------
  
  #pragma omp parallel for
  for(int fID=0;fID<mesh->getNbIntFaces();fID++){
    
    int eTL     = mesh->getIntFacesTag(fID,0);                               
    int eIDL    = mesh->getIntFacesTag(fID,1); 
    int eFNumL  = mesh->getIntFacesTag(fID,2);
    int nbFNPL  = (*fs)[eTL]->getNbFNP();
    const DGIntMatrix* indexL = (*fs)[eTL]->getFRefNodesIndex();
    
    int nbFQP              = (*fs)[eTL]->getNbFQP();
    const DGVector* normal = mesh->getIntFacesNormal(fID);

    int eTR     = mesh->getIntFacesTag(fID,3);                               
    int eIDR    = mesh->getIntFacesTag(fID,4); 
    int eFNumR  = mesh->getIntFacesTag(fID,5);
    int nbFNPR  = (*fs)[eTR]->getNbFNP();
    const DGIntMatrix* indexR = (*fs)[eTR]->getFRefNodesIndex();
    
    // --- get solution ---
    
    DGMatrix UL(nbFNPL,1); for(int i=0;i<nbFNPL;i++) UL(i,0) = (*U)[eTL]((*indexL)(eFNumL,i),eIDL);
    DGMatrix UR(nbFNPR,1); for(int i=0;i<nbFNPR;i++) UR(i,0) = (*U)[eTR]((*indexR)(eFNumR,nbFNPR-1-i),eIDR);     // talking nodes
    
    // --- get solution at quadrature points ---
    
    DGMatrix ULQP(nbFQP,1); ULQP.blas_dgemm(*(*fs)[eTL]->getPhiFQP(),UL,1,0);
    DGMatrix URQP(nbFQP,1); URQP.blas_dgemm(*(*fs)[eTR]->getPhiFQP(),UR,1,0);
        
    // --- copy data to containers ---
    
    DGVector(UFQP[0].touchColumn(fID)).memcpy(DGVector(ULQP.seeColumn(0)));
    DGVector(UFQP[1].touchColumn(fID)).memcpy(DGVector(URQP.seeColumn(0)));
    DGVector(UFNP[0].touchColumn(fID)).memcpy(DGVector(UL.seeColumn(0)));
    DGVector(UFNP[1].touchColumn(fID)).memcpy(DGVector(UR.seeColumn(0)));
      
  }
  
  // --- get solution at boundaries interfaces quadrature points ---------------
  
  #pragma omp parallel for
  for(int bID=0;bID<mesh->getNbBFaces();bID++){
    
    int eTL     = mesh->getBFacesTag(bID,0);                               
    int eIDL    = mesh->getBFacesTag(bID,1); 
    int eFNumL  = mesh->getBFacesTag(bID,2);
    int nbFNPL  = (*fs)[eTL]->getNbFNP();
    const DGIntMatrix* indexL = (*fs)[eTL]->getFRefNodesIndex();
    
    int nbFQP              = (*fs)[eTL]->getNbFQP();
        
    // --- get solution ---
    
    DGMatrix UL(nbFNPL,1); for(int i=0;i<nbFNPL;i++) UL(i,0) = (*U)[eTL]((*indexL)(eFNumL,i),eIDL);
    DGMatrix UR(nbFNPL,1); for(int i=0;i<nbFNPL;i++) UR(i,0) = flux->getBCValue(bID,nbFNPL-1-i);    //nbFNP-1-i to get correct nodes         
    
    // --- get solution at quadrature points ---
    
    DGMatrix ULQP(nbFQP,1); ULQP.blas_dgemm(*(*fs)[eTL]->getPhiFQP(),UL,1,0);
    DGMatrix URQP(nbFQP,1); URQP.blas_dgemm(*(*fs)[eTL]->getPhiFQP(),UR,1,0);
    
    // --- copy data to containers ---
    
    DGVector(UBQP[0].touchColumn(bID)).memcpy(DGVector(ULQP.seeColumn(0)));
    DGVector(UBQP[1].touchColumn(bID)).memcpy(DGVector(URQP.seeColumn(0)));
    
  }
  
}
//------------------------------------------------------------------------------
void Model::computeFUQP(){
  
  // --- get flux/dFlux at elements quadrature points --------------------------
  
  int eT = 0;
  for(SolContainer::const_iterator eTIt=UEQP.begin();eTIt!=UEQP.end();++eTIt,eT++){
    for(int dim=0;dim<2;dim++){
      flux->getFlux(dim,&UEQP[eT],&FUEQP[eT][dim],EQP,eT,0);
      flux->getFluxDerivative(dim,&UEQP[eT],&dFUEQP[eT][dim],EQP,eT,0);
    }
  }
  
  // --- get flux/dFlux at interfaces quadrature points ------------------------
  
  #pragma omp parallel for
  for(int fID=0;fID<mesh->getNbIntFaces();fID++){
    
    int eTL     = mesh->getIntFacesTag(fID,0);                               
    int eIDL    = mesh->getIntFacesTag(fID,1);                               
    int eTR     = mesh->getIntFacesTag(fID,3); 
    int eIDR    = mesh->getIntFacesTag(fID,4);                               
    int nbFQP   = (*fs)[eTL]->getNbFQP();
    int nbFNP   = (*fs)[eTL]->getNbFNP();
    
    for(int dim=0;dim<2;dim++){
      flux->getFlux(dim,&DGMatrix(UFQP[0].seeSubmatrix(0,nbFQP,fID,1)),&DGMatrix(FUFQP[0][dim].touchSubmatrix(0,nbFQP,fID,1)),FQP,eTL,fID);
      flux->getFlux(dim,&DGMatrix(UFQP[1].seeSubmatrix(0,nbFQP,fID,1)),&DGMatrix(FUFQP[1][dim].touchSubmatrix(0,nbFQP,fID,1)),FQP,eTR,fID);
      flux->getFlux(dim,&DGMatrix(UFNP[0].seeSubmatrix(0,nbFNP,fID,1)),&DGMatrix(FUFNP[0][dim].touchSubmatrix(0,nbFNP,fID,1)),FNP,eTL,fID);
      flux->getFlux(dim,&DGMatrix(UFNP[1].seeSubmatrix(0,nbFNP,fID,1)),&DGMatrix(FUFNP[1][dim].touchSubmatrix(0,nbFNP,fID,1)),FNP,eTR,fID);
      flux->getFluxDerivative(dim,&DGMatrix(UFQP[0].seeSubmatrix(0,nbFQP,fID,1)),&DGMatrix(dFUFQP[0][dim].touchSubmatrix(0,nbFQP,fID,1)),FQP,eTL,fID);
      flux->getFluxDerivative(dim,&DGMatrix(UFQP[1].seeSubmatrix(0,nbFQP,fID,1)),&DGMatrix(dFUFQP[1][dim].touchSubmatrix(0,nbFQP,fID,1)),FQP,eTR,fID);
    } 
  }
  
  // --- get flux/dFlux at boundaries interfaces quadrature points -------------
  
  #pragma omp parallel for
  for(int bID=0;bID<mesh->getNbBFaces();bID++){                                 // NOT YET IMPLEMENTED Dirichlet BC
    
    int eTL     = mesh->getIntFacesTag(bID,0);                               
    int eIDL    = mesh->getIntFacesTag(bID,1);                                                              
    int nbFQP   = (*fs)[eTL]->getNbFQP();
    
    for(int dim=0;dim<2;dim++){
      flux->getFlux(dim,&DGMatrix(UBQP[0].seeSubmatrix(0,nbFQP,bID,1)),&DGMatrix(FUBQP[0][dim].touchSubmatrix(0,nbFQP,bID,1)),BQP,eTL,bID);
      flux->getFlux(dim,&DGMatrix(UBQP[1].seeSubmatrix(0,nbFQP,bID,1)),&DGMatrix(FUBQP[1][dim].touchSubmatrix(0,nbFQP,bID,1)),BQP,eTL,bID);
      flux->getFluxDerivative(dim,&DGMatrix(UBQP[0].seeSubmatrix(0,nbFQP,bID,1)),&DGMatrix(dFUBQP[0][dim].touchSubmatrix(0,nbFQP,bID,1)),BQP,eTL,bID);
      flux->getFluxDerivative(dim,&DGMatrix(UBQP[1].seeSubmatrix(0,nbFQP,bID,1)),&DGMatrix(dFUBQP[1][dim].touchSubmatrix(0,nbFQP,bID,1)),BQP,eTL,bID);
    } 
  }
  
}
//------------------------------------------------------------------------------
void Model::computeVolumeTerm(SolContainer* Res){
    
  int eT=0;
  for(SolContainer::const_iterator eTIt=Res->begin();eTIt!=Res->end();++eTIt,eT++){
        
    int nbEQP  = (*fs)[eT]->getNbEQP();
    int nbElmt = mesh->getNbElmt(eT);

    SolContainer FluxUVQP(2,DGMatrix(nbEQP,nbElmt)); 
    
    // --- Loop over the elements ---

    //#pragma omp parallel for
    for(int eID=0;eID<nbElmt;eID++){
      
      DGMatrix* invJ = mesh->getInvJ(eT,eID);
    
      // --- compute fuv ---
      
      for(int dimUV=0;dimUV<2;dimUV++){
        for(int dimXY=0;dimXY<2;dimXY++){
          DGVector(FluxUVQP[dimUV].touchColumn(eID)).scaleAndAdd((*invJ)(dimXY,dimUV),DGVector(FUEQP[eT][dimXY].seeColumn(eID)));
        }
      }

    }

    // --- redistribute to nodal points ---
    
    (*Res)[eT].set_all(0.);
    for(int dimXY=0;dimXY<2;dimXY++) (*Res)[eT].blas_dgemm_omp_col(*(*fs)[eT]->getFqduv(dimXY),FluxUVQP[dimXY],1,1);
             
  }
  
  
}
//------------------------------------------------------------------------------
void Model::computeInterfaceTerm(SolContainer* Res){
  
  for(int fID=0;fID<mesh->getNbIntFaces();fID++){
    
    int eTL      = mesh->getIntFacesTag(fID,0);                               
    int eIDL     = mesh->getIntFacesTag(fID,1); 
    int eFNumL   = mesh->getIntFacesTag(fID,2);
    int nbFNPL   = (*fs)[eTL]->getNbFNP();
    double detJL = mesh->getDetJ(eTL,eIDL);
    const DGIntMatrix* indexL = (*fs)[eTL]->getFRefNodesIndex();
    
    int nbFQP              = (*fs)[eTL]->getNbFQP();
    int nbFNP              = (*fs)[eTL]->getNbFNP();
    const DGVector* normal = mesh->getIntFacesNormal(fID);
    
    int eTR      = mesh->getIntFacesTag(fID,3);                               
    int eIDR     = mesh->getIntFacesTag(fID,4); 
    int eFNumR   = mesh->getIntFacesTag(fID,5);
    int nbFNPR   = (*fs)[eTR]->getNbFNP();
    double detJR = mesh->getDetJ(eTR,eIDR);
    const DGIntMatrix* indexR = (*fs)[eTR]->getFRefNodesIndex();
    
    // --- get numerical flux at interface quadrature points ------------------- 
    
    DGMatrix FluxQPn(nbFQP,1);
    flux->getNumericalFlux(&DGMatrix(UFQP[0].seeSubmatrix(0,nbFQP,fID,1)),
                           &FUFQP[0],
                           &dFUFQP[0], 
                           &DGMatrix(UFQP[1].seeSubmatrix(0,nbFQP,fID,1)),
                           &FUFQP[1],
                           &dFUFQP[1],
                           &FluxQPn,
                           normal,
                           fID);
    
    // --- redistribute residu to nodes ----------------------------------------
    
    DGMatrix FluxNPn(nbFNP,1);
    FluxNPn.blas_dgemm(*(*fs)[eTL]->getFquv(),FluxQPn,1,0);
    
    for(int i=0;i<nbFNP;i++){
      (*Res)[eTL]((*indexL)(eFNumL,i),eIDL) -= 0.5*(*normal)(2)/detJL*FluxNPn(i,0);
      (*Res)[eTR]((*indexR)(eFNumR,i),eIDR) += 0.5*(*normal)(2)/detJR*FluxNPn(nbFNP-1-i,0);  //nbFNP-1-i to get correct nodes
    }

  }

}
//------------------------------------------------------------------------------
void Model::computeBoundaryTerm(SolContainer* Res){
  
  #pragma omp parallel for
  for(int bID=0;bID<mesh->getNbBFaces();bID++){
    
    int eTL                   = mesh->getBFacesTag(bID,0);                               
    int eIDL                  = mesh->getBFacesTag(bID,1); 
    int eFNumL                = mesh->getBFacesTag(bID,2);
    double detJL              = mesh->getDetJ(eTL,eIDL);
    int nbFQP                 = (*fs)[eTL]->getNbFQP();
    int nbFNP                 = (*fs)[eTL]->getNbFNP();
    const DGVector* normal    = mesh->getBFacesNormal(bID);
    const DGIntMatrix* indexL = (*fs)[eTL]->getFRefNodesIndex();

    // --- get numerical flux at interface quadrature points ------------------- 
    
    DGMatrix FluxQPn(nbFQP,1);
    flux->getNumericalFlux(&DGMatrix(UBQP[0].seeSubmatrix(0,nbFQP,bID,1)),
                           &FUBQP[0],
                           &dFUBQP[0], 
                           &DGMatrix(UBQP[1].seeSubmatrix(0,nbFQP,bID,1)),
                           &FUBQP[1],
                           &dFUBQP[1],
                           &FluxQPn,
                           normal,
                           bID);
    
    // --- redistribute residu to nodes ----------------------------------------
    
    DGMatrix FluxNPn(nbFNP,1);
    FluxNPn.blas_dgemm(*(*fs)[eTL]->getFquv(),FluxQPn,1,0);
    
    for(int i=0;i<nbFNP;i++) (*Res)[eTL]((*indexL)(eFNumL,i),eIDL) -= 0.5*(*normal)(2)/detJL*FluxNPn(i,0);
    
  }
  
}
//------------------------------------------------------------------------------
void Model::setBilinearForm(){
  
  //ATTENTION: hard coded for p1!
  
  //has to be multiply by elmt area
  bE[0](0,0) =  1./2.; bE[0](0,1) = -1./4.; bE[0](0,2) = -1./4.;
  bE[0](1,0) = -1./4.; bE[0](1,1) =  1./2.; bE[0](1,2) = -1./4.;
  bE[0](2,0) = -1./4.; bE[0](2,1) = -1./4.; bE[0](2,2) =  1./2.;
  
  //has to be multiply by face length
  bF(0,0) =  1.   ; bF(0,1) =  0.   ; bF(0,2) = -2./3.; bF(0,3) = -1./3.;
  bF(1,0) =  0.   ; bF(1,1) =  1.   ; bF(1,2) = -1./3.; bF(1,3) = -2./3.;
  bF(2,0) = -2./3.; bF(2,1) = -1./3.; bF(2,2) =  1.   ; bF(2,3) =  0.   ;
  bF(3,0) = -1./3.; bF(3,1) = -2./3.; bF(3,2) =  0.   ; bF(3,3) =  1.   ;
 
}
//------------------------------------------------------------------------------
void Model::computeArtViscousVolumeTerm(const SolContainer* U,SolContainer* Res){
  
  int eT=0;
  for(SolContainer::const_iterator eTIt=Res->begin();eTIt!=Res->end();++eTIt,eT++){
    
    int nbEQP  = (*fs)[eT]->getNbEQP();
    int nbENP  = (*fs)[eT]->getNbENP();
    
    // --- Loop over the elements --- ResE_i = Sum_j -nuE*u_j*bE(j,i)
    
    #pragma omp parallel for
    for(int eID=0;eID<mesh->getNbElmt(eT);eID++){
      for(int i=0;i<nbENP;i++){
        DGMatrix((*Res)[eT].touchSubmatrix(i,1,eID,1)).blas_dgemm(DGMatrix(bE[eT].seeSubmatrix(i,1,0,nbENP)),
                                                                  DGMatrix((*U)[eT].seeSubmatrix(0,nbENP,eID,1)),
                                                                  -nuE[eT](eID),1);
                                                                  //*1 because |K|=detJ is not included into bE, but detJ is killed by detJ from invMass!
      }
    }
    
  }
  
}
//------------------------------------------------------------------------------
void Model::computeArtViscousInterfaceTerm(const SolContainer* U,SolContainer* Res){
  
  // --- Loop over the interfaces --- ResF_i = -nuF*u_i*bF(i,i) + Sum_j_ext -nuF*u_j*bF(j,i)

  for(int fID=0;fID<mesh->getNbIntFaces();fID++){
    
    int eTL      = mesh->getIntFacesTag(fID,0);                               
    int eIDL     = mesh->getIntFacesTag(fID,1); 
    int eFNumL   = mesh->getIntFacesTag(fID,2);
    double detJL = mesh->getDetJ(eTL,eIDL);
    const DGIntMatrix* indexL = (*fs)[eTL]->getFRefNodesIndex();
    
    int nbFQP              = (*fs)[eTL]->getNbFQP();
    int nbFNP              = (*fs)[eTL]->getNbFNP();
    const DGVector* normal = mesh->getIntFacesNormal(fID);
    
    int eTR      = mesh->getIntFacesTag(fID,3);                               
    int eIDR     = mesh->getIntFacesTag(fID,4); 
    int eFNumR   = mesh->getIntFacesTag(fID,5);
    double detJR = mesh->getDetJ(eTR,eIDR);
    const DGIntMatrix* indexR = (*fs)[eTR]->getFRefNodesIndex();
        
    // --- get solution at face nodal points ---
    
    DGMatrix UL(nbFNP,1);
    DGMatrix UR(nbFNP,1);
    
    for(int i=0;i<nbFNP;i++){
      UL(i,0) = (*U)[eTL]((*indexL)(eFNumL,i),eIDL);
      UR(i,0) = (*U)[eTR]((*indexR)(eFNumR,nbFNP-1-i),eIDR);                    // nbFNP-1-i to have talking nodes
    }
        
    // --- compute interface residual term ---                                    
    
    DGMatrix tmpFResL(nbFNP,1);
    DGMatrix tmpFResR(nbFNP,1);

    for(int i=0;i<nbFNP;i++){ 							// multiply by normal because |F|=normal(2) and is not included into bF
      tmpFResL(i,0) = -UL(i,0)*nuF(fID)*(*normal)(2)*bF(i,i);
      tmpFResR(i,0) = -UR(i,0)*nuF(fID)*(*normal)(2)*bF(i,i);
      for(int j=0;j<nbFNP;j++){
        tmpFResL(i,0) += -UR(j,0)*nuF(fID)*(*normal)(2)*bF(nbFNP+j,i);
        tmpFResR(i,0) += -UL(j,0)*nuF(fID)*(*normal)(2)*bF(nbFNP+j,i);        
      }
    }
    
    // --- redistribute to nodes ---
    
    for(int i=0;i<nbFNP;i++){
      (*Res)[eTL]((*indexL)(eFNumL,i),eIDL) += 1./detJL*tmpFResL(i,0);             // 1/detJ from the invMass
      (*Res)[eTR]((*indexR)(eFNumR,i),eIDR) += 1./detJR*tmpFResR(nbFNP-1-i,0);     // nbFNP-1-i to get back to elmt ordering
    }
    
  }
  
}
//------------------------------------------------------------------------------
void Model::computeViscosity(const SolContainer* U){
    
  int nbF = mesh->getNbIntFaces();
  SolContainer tmpIFTerm;
  
  // -- Compute nuE ------------------------------------------------------------

  int eT=0;
  for(SolContainer::const_iterator eTIt=U->begin();eTIt!=U->end();++eTIt,eT++){
    
    int nbEQP  = (*fs)[eT]->getNbEQP();
    int nbENP  = (*fs)[eT]->getNbENP();
    int nbFNP  = (*fs)[eT]->getNbFNP();
    int nbElmt = mesh->getNbElmt(eT);
    DGMatrix tmpERes(nbENP,nbENP);
    for(int LR=0;LR<2;LR++) tmpIFTerm.push_back(DGMatrix(nbFNP,nbF*nbFNP));    
    
    #pragma omp parallel for
    for(int eID=0;eID<nbElmt;eID++){
      
      SolContainer dFluxUVQP;
      for(int dim=0;dim<2;dim++) dFluxUVQP.push_back(DGMatrix(nbEQP,1));
        
      DGMatrix* invJ = mesh->getInvJ(eT,eID);
      double    detJ = mesh->getDetJ(eT,eID);
      
      // --- compute fuv_prime ---
      
      for(int dimUV=0;dimUV<2;dimUV++){
        for(int dimXY=0;dimXY<2;dimXY++){
          DGVector(dFluxUVQP[dimUV].touchColumn(0)).scaleAndAdd((*invJ)(dimXY,dimUV),DGVector(dFUEQP[eT][dimXY].seeColumn(eID)));
        }
      }
      
      // --- compute vol term: IntE (f'.dphi_j)*phi_i ---
      
      tmpERes.set_all(0.);                                                      // tmpERes = (j,i)
      for(int dimXY=0;dimXY<2;dimXY++){
        for(int i=0;i<nbENP;i++){
          DGMatrix(tmpERes.touchSubmatrix(0,nbENP,i,1)).blas_dgemm(*(*fs)[eT]->getdPhiUVPhiEQP(dimXY,i),
                                                                   dFluxUVQP[dimXY],
                                                                   detJ,1);
        }
      }
      
      // --- compute interface term: IntF (S.n)*phi_j*phi_i --- S = (numf-f_int)/jump 

      for(int f=0;f<mesh->getNbElmtFaces(eT);f++){
        
        int fID      = mesh->getElmtFaceID(eT,eID,f);
        int LorR     = mesh->getElmtFaceID(eT,eID,f,1);
        if(LorR==3) continue;                                                   //ATTENTION: BC is not yet included
        
        int L        = LorR;
        int eTL      = eT;                               
        int eIDL     = eID; 
        int eFNumL   = f;
        int nbFNPL   = (*fs)[eTL]->getNbFNP();
        double detJL = detJ;
        const DGIntMatrix* indexL = (*fs)[eTL]->getFRefNodesIndex();
        
        int nbFQP              = (*fs)[eTL]->getNbFQP();
        int nbFNP              = (*fs)[eTL]->getNbFNP();
        DGVector normal(3); normal.memcpy(*mesh->getIntFacesNormal(fID));       // normal sign is function of LorR!!!
        if(LorR==1){normal(0) = -normal(0);normal(1) = -normal(1);}
          
        int R        ;LorR==0? R=1:R=0;
        int eTR      = mesh->getIntFacesTag(fID,3-3*LorR);                               
        int eIDR     = mesh->getIntFacesTag(fID,4-3*LorR); 
        int eFNumR   = mesh->getIntFacesTag(fID,5-3*LorR);
        int nbFNPR   = (*fs)[eTR]->getNbFNP();
        double detJR = mesh->getDetJ(eTR,eIDR);
        const DGIntMatrix* indexR = (*fs)[eTR]->getFRefNodesIndex();

	// --- get jump solution at quadrature points ---                       // ordering is function of LorR, so take care when you refer to elmt node!
        
        DGMatrix jumpQP(nbFQP,1);
        if(LorR==0){
          for(int q=0;q<nbFQP;q++) jumpQP(q,0) = UFQP[L](q,fID)-UFQP[R](q,fID);
        }else{
          for(int q=0;q<nbFQP;q++) jumpQP(q,0) = UFQP[L](nbFQP-1-q,fID)-UFQP[R](nbFQP-1-q,fID);
        }
        
        // --- get numerical flux at quadrature points ---
        
        DGMatrix FluxQPn(nbFQP,1);
        flux->getNumericalFlux(&DGMatrix(UFQP[L].seeSubmatrix(0,nbFQP,fID,1)),
                               &FUFQP[L],
                               &dFUFQP[L], 
                               &DGMatrix(UFQP[R].seeSubmatrix(0,nbFQP,fID,1)),
                               &FUFQP[R],
                               &dFUFQP[R],
                               &FluxQPn,
                               &normal,
                               fID);
        
        if(LorR==1){
          DGVector tmp(nbFQP);
          tmp.memcpy(DGVector(FluxQPn.touchColumn(0)));
          for(int q=0;q<nbFQP;q++) FluxQPn(nbFQP-1-q,0) = tmp(q);
        }
        
        // --- get flux ----------
        
        DGMatrix tmpFRes(nbFNP,nbFNP);
        
        for(int i=0;i<nbFNP;i++){
                                                                                // S_QP_ni = (numF_QPn - F_QPn)/jump_QP  where F_QPn is the interior flux (L)
          DGMatrix S(nbFQP,1);          
          for(int q=0;q<nbFQP;q++){
            if(fabs(jumpQP(q,0))<1e-10) {S(q,0) = 0.; continue;}
            if(LorR==0)  {S(q,0) = (FluxQPn(q,0)-(normal(0)*FUFQP[L][0](q,fID)+normal(1)*FUFQP[L][1](q,fID)))/jumpQP(q,0);}
            else         {S(q,0) = (FluxQPn(q,0)-(normal(0)*FUFQP[L][0](nbFQP-1-q,fID)+normal(1)*FUFQP[L][1](nbFQP-1-q,fID)))/jumpQP(q,0);}
          }
          
          // --- add interface term: IntF Sn*phi_j*phi_i ---
          
          DGMatrix(tmpFRes.touchSubmatrix(0,nbFNP,i,1)).blas_dgemm(*(*fs)[eT]->getPhiPhiFQP(i),S,0.5*normal(2),0);		//*0.5*normal because detF = 0.5*n

          // --- redistribute to Face nodes of elmt eIDL ---
          
          for(int j=0;j<nbFNP;j++) tmpERes((*indexL)(eFNumL,j),(*indexL)(eFNumL,i)) += tmpFRes(j,i);
        
        }
        
        // --- save interface term ---
        
        DGMatrix(tmpIFTerm[L].touchSubmatrix(0,nbFNP,fID*nbFNP,nbFNP)).memcpy(tmpFRes);
        
      }
      
      // --- get nuE --------
      
      nuE[eT](eID) = 0.;									//reinitiliaze nuE
      
      for(int i=0;i<nbENP;i++){
        for(int j=0;j<nbENP;j++){
          nuE[eT](eID) = max(nuE[eT](eID),max(tmpERes(j,i),0.)/fabs(detJ*bE[eT](i,j)));      // detJ because area is not included into bE
        }
      }      
      
            
    }
    
  }
  
  
  // -- compute nuF ------------------------------------------------------------
  
  #pragma omp parallel for  
  for(int fID=0;fID<mesh->getNbIntFaces();fID++){
    
    int eTL   = mesh->getIntFacesTag(fID,0);                               
    int nbFNP = (*fs)[eTL]->getNbFNP();
    const DGVector* normal = mesh->getIntFacesNormal(fID);
    
    // --- get solution at face nodal points ---
    int a = 0;
    int b = 0;
    nuF(fID) = 0.; 										//reinitiliaze nuF
      
    for(int k=0;k<2;k++){
      if(k==0) {a=0;b=1;}
      else     {a=1;b=0;}
      for(int i=0;i<nbFNP;i++){
        for(int j=0;j<nbFNP;j++){
          nuF(fID) = max(nuF(fID),max(0.,-tmpIFTerm[k](j,fID*nbFNP+i))/fabs((*normal)(2)*bF(a*nbFNP+j,b*nbFNP+i)));
        }
      }
    }
    
  }

}

